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It is hard to solve ill-posed problems, as calculated temperatures are very sensitive to errors made while calculat- 


ing “measured” temperatures or performing real-time measurements. The errors can create temperature oscillation, 


which can be the cause of an unstable solution. In order to overcome such difficulties, a variety of techniques 


have been proposed in literature, including regularization, future time steps and smoothing digital filters. In this 


paper, the Tikhonov regularization is applied to stabilize the solution of the inverse heat conduction problem. The 


impact on the inverse solution stability and accuracy is demonstrated. 
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Introduction 


The thermal analysis is necessary for the design of 
components operating in high temperatures. One of the 
difficulties in the thermal analysis is the definition of all 
boundary conditions, especially the convection boundary 
condition. The problem can be solved numerically, ana- 
lysing the flow phenomena occurring when the fluid 
comes into contact with individual components [1]. An 
alternative method of the temperature distribution deter- 
mination is the search for the solution to the inverse heat 
conduction problem in the device under analysis. 

Inverse problems, in contrast to direct problems, are 
characterized by an unknown part of boundary conditions 
[2-4, 9], some undefined material properties [5] or inde- 
terminate heat sources [6]. They happen in engineering, 
because defining some boundary conditions, material 
properties or heat sources is often tough. These problems 
can be solved by performing additional temperature 
measurements instead of the unknown boundary condi- 
tion or the parameters mentioned above. 
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Formulation of the method 


Using the semi-discrete control volume method to de- 
termine the temperature distribution in a plate, as presented 
in Fig. 1, a system of non-linear differential equations can 
be obtained. 
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Fig. 1 Plate with thickness L with an unknown and known 
boundary condition 
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Nomenclature 
c specific heat, J/(kgK) T temperature, °C or K 
k thermal conductivity, W/mK) q heat flux, W/m? 
t time, s 


Fig. 2 Staircase approximation of temperature or heat flux at 
the plate front surface 


If the temperature and the heat flux at the body surface 
are approximated by a continuous-time staircase function 
both at the plate front (cf. Fig. 2) and back surface 

u(t) =u(k) 

w(t) = w(k) 

the following system of equations describes temperature 
distribution 


kAt<t<(k+1)At,k€(0,..N) (1) 


T=AT+Bu+ Ww (2) 
with the initial conditions 
T, (to) 
T, (to) 
Toa 6) 
Ty, (to) 


Ng 

where T is the Ng dimensional vector, which describes 
temperature transients at Ng nodes of a discrete space 
model (state vector), time variable t varies in the range of 
t, <t<t,,,, A is the square matrix of Ng x Ng dimen- 


sions, B and W are vectors of the Ng dimensions, u(t) is 
the unknown boundary condition, w(t) is the known 
boundary condition. Both boundary conditions can be 
assumed in the form of the first or the second kind 
boundary condition. The solution of (2) can be written in 
a matrix form [7] 


T(t)= ep t) + call e *“Bu(s)ds 


i (4) 
+e f, e ““ww/(s)ds 
0 


Replacing t with to+At in equation (4), the following is 
obtained: 


T(t) + At) =e*“T (t,) 


+ ellotACJA tat -sA Bu (s)ds (5) 


to 
otto Ata oe Aww (s)ds 
0 
As noted in (1), u(s) and w(s) are held constant 
throughout the interval at the value of u(kAt) and w(kAt), 
respectively. Therefore, equation (5) may be rewritten as 
T(t) + At) = F(At)T(t))+G(At)u(t))+H(At)w(t)) (6) 
Where 


+ 


F (At) = e^“ (7) 
G(At) = | e^dsB (8) 
H(At) = i. eSdsw (9) 


Matrix G(4t) is given by 


G(At)=G= feo gds 
° (10) 


0 A 
=-| e^ Bdr =| ‘eATdTB 
At 0 


where 
T =ty+At—s . Matrix H can be calculated in a similar 
way. Expanding e* into the McLaurin series in the same 
way as function e*, F(At)= F can be written as 
2 3 
AAt AAt 
(aar) + nM +... (11) 
2! 3! 
Introducing e^™ into (8) and integrating (8) gives 
G(At)= G 


F=I+AA^At+ 


2 3 
G= GAN py s pae (12) 


2! 3! 4! 


Matrix H(At)= H can be calculated in a similar way. 


2 3 
galg, Ade, (AAO? (AAt) su [war (13) 


2! 3! 4! 


Using equation (6) and omitting At, the temperature in 
time t; is obtained 


T(t,)=T, =FT, + Gu, + Hwy (14) 
Similarly, for time point 2 
T, =FT, + Gu, + Hw, (15) 


Substituting (14) into (15) gives 
T, = F[FT) + Gu) +Hw)|]+Gu,+Hw, (16) 
after transformation 
T, = FT + FGu, + FHw, + Gu, + Hw, (17) 
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and generally for point i 
T; = FİT, + F} ‘Gu, + F} Hw, + F'°Gu, 
+F} Hw; +...+FGu,_, +FHw,_ (18) 
+ Gu,_; + Hw;_;, 


T, =F'T,+ F Gu +$ F Hw, (19) 
j=0 j=0 
In order to determine the temperature at a selected 
spatial point L, the following row vector is introduced 


V; = (OseosOp oo) y, (20) 


where a,=1. The temperature at spatial point L and at 
time i can be written as 


i-l i-l 
T i =v|F'T +X F Gu +Y F Hw CD 
j=0 j=0 
Let us assume that temperature is measured at nodes: 
Li, La, p Ly,- The aim of the proposed inverse method 
is to calculate the unknown boundary condition u(t) so 
that the calculated temperatures T agree, within certain 
limits, with experimentally measured temperatures f. 


fii =T July, i=1,..,(N+1) (22) 
If there is just one measurement point at the insulated 
surface, L=L4. 


fi, SVa [FT + Guy + Hwg] 
fu2= Vu (FT, + FGu, + Gu, +FHw, + Hw] 
EO Rm (23) 


N 
N+ N-j 
faa = Vn r "Ty + J F" Gu; 
j=0 


N 
+) et 
j=0 


The above system of equations (23) can be written in a 
matrix form 


f-c=ZU (24) 
where 
fra 
i- oo (25) 
fr N+ 


Va [FT) + Hw | 
c= a (26) 


N . 
Va | FT, +>) F* /Hw, 
j=0 
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Up 
Uy, =4 (27) 
Uy 
VG 0 . 0 
V FG VG 0 
Z= P (28) 


v,FXG v FG ... vG 


1 


The least squares method is used to satisfy the condi- 
tion expressed in (22) 
Nr N+ 


> > (Tan fu) =min (29) 
j=l i=l 
This can be written in the following matrix form 
S =(ZU-f+c)' (ZU-f+c)=min (30) 
The solution which satisfies condition (30) is 
u=(z"z) 2" (fc) (31) 


The determinant of matrix Z is close to zero and cal- 
culations according to (31) could often lead to an 
ill-conditioned algorithm. Hence, it is necessary to inves- 
tigate the problem, estimate its range and look for ways 
of reducing it. 

One of the possibilities of improving the inverse 
problem solution stability is the Tikhonov regularization 
[8]. Tikhonov proposes parallel minimization of the qua- 
dratic measure of the solution vector U. 

S =|ZU-f + c|? +2 Ul? = min (32) 

Condition (32) leads to 

U=(z"Z+Al) ‘2! (fc) (33) 

where I is a unit matrix. Here, the mechanism of im- 
proving the conditioning of matrix Z'Z can be observed. It 
consists in adding constants à to the matrix diagonal. The 
weighting coefficient à should be relatively small so as not 
to change the main purpose of equation (30) — minimiza- 
tion of the square residuals — but it cannot be too small, 
especially in the case of strongly ill-conditioned equation 
(31). The optimal choice of à will be discussed in the next 
section. 

Total errors in the identified boundary condition u(t) can 
be estimated by comparing calculated values u, with exact 


values u,;” for i=0, N using the following formula: 


(34) 


Similarly, the total error in the calculated temperature 
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distribution can be obtained 


Ng N+ 2 
ZAT -T#) 


{== 35 


Analytical Solution 


Let us assume that an infinite plate with thickness L, as 
presented in Fig. 1, is heated by a heat flux on its front sur- 
face and that the plate back surface is insulated. The initial 
temperature of the plate is constant and equals Tp. Material 
properties are assumed as temperature-independent. The 
temperature distribution in the plate is described by the heat 
conduction equation 


2 
a =a = O<x<L (36) 
ot Ox 
the boundary conditions 
tı<t<t 
_ôT| Jim i 2 (37) 
OX |,-9 (0 t<t; or t2t 
aT) 26 (38) 
Ox Xx=L 
and the initial condition 
T(x,0)=0 O<x<L (39) 


The plate temperature distribution can be calculated us- 
ing the superposition principle [4]. 
T(x,t)=0 O<t<ti (40) 
T (x,t) =9(x,t-t) ti<t<h (41) 
T (x,t)=9(x,t-t)-9(x,t-t) t>t (42) 
Where 


2, 
son eE 4 2 es 


k |e 2 L 3 


2 %1 t 
— 2Y eos[ne *) exp (=z z) 


T 


Numerical simulations 


The plate presented in Fig. 1 with thickness L=0.1 m 
is insulated on its back surface and is subjected to the 
boundary condition of the second kind (37) with para- 
meters: t1=730 s, t= 1230 s, qm=350 kW/m’. 

The following material properties are adopted for the 
solution: k=40 W/(mK), co=4.0E+6 J/(m°K). 

The inverse problem is then formulated without speci- 
fying the boundary condition on the plate front surface 
with the measured temperature transient on the back insu- 
lated surface. The heat flux transient at the front surface 
of the plate with thickness L=0.1 m and the temperature 
distribution on the thickness are determined based on the 
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measured temperature history at the plate back surface, 
which is insulated. “Exact measurement data” at the back 
surface are obtained from the analytical solution (40-43). 
The plate thickness is divided into 6 control volumes 
(N:=6) and the time step At=S0 s is used. 

The algorithm for calculating the unknown boundary 
condition vector U, which is based on equation (31), is 
unstable. The results are presented in Fig. 3. 

The solution formulated according to equation (33) 
using the Tikhonov regularization method makes it pos- 
sible to find the unknown boundary condition. The 
smallest error is found with the weighting coefficient A= 
5.7E-10. The calculated heat flux at the plate front sur- 
face and the determined temperature transients are 
presented in Fig. 4 and Fig. 5. Total errors in the identified 
boundary condition and in the temperature distribution 
equal to ©,=4.49E+04 W/m? and £= 4.41°C, respectively. 
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Fig. 3 Estimated unknown boundary condition at the plate 
front surface according to (31) without regularization 
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Fig. 4 Estimated unknown boundary condition at the plate 
front surface according to (33) using the regularization 
method 
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Fig. 5 Comparison between the exact temperature histories 


and those calculated by means of the inverse method 
using regularization 


In order to make the test more realistic, “noisy meas- 
ured data” are obtained by adding normal random errors 
+10°C of zero mean to “exact measured data”. The small- 


est error is found with the weighting coefficient A=1.7E-9. 


It is an order of magnitude larger than for undisturbed 
data. The calculated heat flux at the plate front surface 
and the determined temperature transients are presented 
in Fig. 6 and Fig. 7. Total errors in the identified boundary 
condition and in the temperature distribution equal X= 
4.86E+04 W/m? and X= 7.04°C, respectively. Despite 
large measurement errors, the total error in the identified 
boundary condition is comparable with the identification 
based on undisturbed data. 
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Fig. 6 Estimated unknown boundary condition at the plate 
front surface according to (33) using the regularization 


method based on “noisy measured data” 
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Fig. 7 Comparison between the exact temperature histories 
and those calculated by means of the inverse method 


using regularization based on “noisy measured data” 
Conclusions 


This paper presents a method of solving the inverse heat 
conduction problem. The problem is ill-posed and its solu- 
tion is unstable. In order to improve the solution stability, 
the Tikhonov regularization is proposed. 

The proposed method is applied to identify the heat flux 
at the front surface of a thick plate based on the measured 
temperature history at the plate back surface, which is 
insulated. The “exact measurement data” at the back sur- 
face are obtained from the presented analytical solution. 
In order to make the test more realistic, “noisy measured 
data? are obtained by adding normal random errors 
+10°C of zero mean to the “exact measured data”. 

The paper also presents the mechanism of improving the 
equation system conditioning. The optimal values of the 
weighting coefficient à are determined both for the "exact" 
and the “noisy measured data". It is an order of magni- 
tude larger for disturbed data than for undisturbed data. 
The numerical verification of the presented inverse method 
applied together with the Tikhonov regularization demon- 
strates its good accuracy and stability. 
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